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ABSTRACT 

This  report  investigates  an  approach  to  characterising  Ground  Probing  Radar  (GPR) 
backscatter  echoes  from  landmines  using  linear  combinations  of  exponentially  damped 
sinusoids.  The  GPR  signatures  of  surrogate  landmines  and  PVC  cylinders  buried  in 
dry  sand  are  measured  using  an  impulse  radar  system  with  centre  frequency  of  1.4  GHz 
and  a  90%  bandwidth.  The  GPR  signal  parameters  are  represented  as  sets  of  complex 
poles  computed  from  a  series  of  neighbouring  signatures  recorded  over  each  target  type. 
The  algorithm  proposed  by  Kumaresan  and  Tufts  which  uses  backward  linear  prediction 
and  the  low-rank  data  matrix  approximation  based  on  singular  value  decomposition  is 
applied  to  this  computation. 

The  performance  of  the  Kumaresan  and  Tufts  (KT)  algorithm  is  compared  with  that 
of  the  Prony  method  when  both  techniques  are  applied  to  modelling  simulated  signals. 
It  is  concluded  that  the  KT  method  provides  more  stable  pole  estimates. 

Two  approaches  to  determining  the  order  of  the  model  are  examined  and  compared 
for  simulated  and  real  data. 

The  results  show  that  the  poles  corresponding  to  different  target  types  form  clusters 
in  the  two-dimensional  a-/  space  (where  a  is  the  pole  damping  factor  and  /  is  the  pole 
frequency).  This  indicates  that  these  pole  clusters  can  be  used  for  the  recognition  of 
landmines. 
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An  Approach  to  Characterising  Ground  Probing  Radar  Target 
Echoes  for  Landmine  Recognition 


EXECUTIVE  SUMMARY 

Modern  anti-personnel  landmines  have  a  minimal  metallic  content  which  makes  detection  using 
electromagnetic  induction  (EMI)  detectors  very  difficult,  if  not  impossible.  To  detect  even  a  low- 
metallic-content  mine  the  detector  threshold  must  be  set  so  low  that  many  false  alarms  occur. 
DSTO  is  investigating  the  use  of  Ground  Probing  Radar  (GPR)  as  a  means  of  detecting  non- 
metallic  landmines  and  identifying  mines  with  a  metal  content.  The  GPR  used  is  a  centimetre 
resolution  short  range  radar.  Target  responses  consist  of  a  set  of  reflections  from  dielectric  surfaces 
in  the  ground  and  in  targets  as  well  as  the  metal  components  within  targets. 

This  report  investigates  an  approach  to  characterising  GPR  backscatter  signals  from  landmines 
as  linear  combinations  of  exponentially  damped  sinusoids.  The  targets  used  in  the  experiments  are 
surrogate  landmines  and  PVC  cylinders  of  various  sizes.  The  targets  are  buried  in  dry  sand  and 
the  measurements  were  made  using  a  1.4  GHz  centre  frequency  antenna  with  90%  bandwidth. 

Two  methods  for  computing  the  signal  parameters  (poles)  are  compared  using  the  simulated 
data  contaminated  by  white  Gaussian  noise:  the  Prony  method  and  a  method  proposed  by  Ku- 
maresan  and  Tufts  (KT).  The  KT  algorithm,  which  applies  backward  linear  prediction  and  the 
low-rank  data  matrix  approximation  based  on  singular  value  decomposition,  is  found  to  give  more 
stable  pole  estimates.  This  method  is,  therefore,  applied  to  model  the  real  target  signatures  where 
the  signal  poles  are  computed  using  a  series  of  neighbouring  GPR  signatures  recorded  over  each 
target  type.  Both  raw  and  background-removed  GPR  target  signatures  were  used. 

Two  different  approaches  for  the  model  order  selection  are  also  investigated:  1)  an  information- 
theoretic  approach  and  2)  an  approach  based  on  the  thresholding  of  singular  values  of  the  data 
matrix.  While  the  information-theoretic  criterion  was  successful  for  the  simulated  data,  the  simpler 
thresholding  approach  has  been  found  to  give  better  results  for  modelling  real  signals.  The  poles 
computed  using  the  KT  algorithm  and  this  latter  approach  were  stable  for  a  range  of  threshold 
values. 

The  results  of  this  initial  study  show  that  the  poles  corresponding  to  different  target  types  form 
clusters  in  the  two-dimensional  a-f  space  (where  a  is  the  pole  damping  factor  and  /  is  the  pole 
frequency).  It  is  therefore  expected  that  such  poles  can  be  used  for  the  recognition  of  landmines. 
However,  further  research  is  needed  to  confirm  these  results. 
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1  Introduction 

Detection  and  identification  of  buried  landmines  is  an  important  yet  difficult  task. 
Standard  countermine  methods  that  use  electromagnetic  induction  (EMI)  can  easily  lo¬ 
cate  metallic  objects.  However,  most  contemporary  landmines  are  characterised  by  a  low 
metallic  content  and  EMI  detectors  often  fail  to  locate  such  devices.  Recent  developments 
in  Ground  Probing  Radar  (GPR)  which  utilise  short  duration  pulses  (0.1-1  ns)  offer  a 
promising  new  approach  to  this  problem.  In  this  situation  antenna  and  target  are  located 
in  half-spaces  having  different  electromagnetic  properties.  The  backscattered  echoes  which 
are  usually  distorted  by  multiple  scattering  are  used  for  target  identification  (see  Daniels 
et  al  [1]  and  Daniels  [2]  for  a  comprehensive  review  of  current  GPR  techniques  and 
their  applications).  A  serious  drawback  of  the  GPR  system  is  that  it  is  limited  in  ability 
to  discriminate  between  landmine  and  non-landmine  echoes.  Advanced  signal  processing 
and  target  classification  techniques  are  sought  with  the  aim  of  solving  this  problem.  It  is 
also  expected  that  improved  detection/recognition  performance  can  be  achieved  through 
combining  different  types  of  sensors,  such  as  EMI  detector,  GPR,  and  infrared  imagery. 

This  report  is  exclusively  concerned  with  the  information  that  can  be  extracted  from 
the  GPR  backscattered  target  echoes.  The  approach  is  to  model  separate  backscattered 
GPR  echoes  from  a  landmine  target  as  a  linear  combination  of  exponential  functions 
exp(st).  The  complex  parameters  (i.e.,  poles)  s  are  expected  to  be  characteristic  for  each 
target-type,  and  invariant  with  respect  to  target  orientation.  This  approach  was  first 
applied  to  the  problem  of  free  space  scattering  from  resonant  bodies  [3].  Chan  et  al.  [4,  5] 
have  demonstrated  that  similar  technique  can  be  successfully  applied  to  identifying  buried 
plastic  (dielectric)  mine-like  targets.  Other  approaches  include  extinction-pulse  (E-pulse) 
or  kill-pulse  (K-pulse)  methods  [6,  7,  8]  applied  to  free-space  scattering,  and  methods 
utilising  artificial  neural  networks  [9,  10]  and  time- frequency  analysis  [11]  as  applied  to 
the  GPR  backscattered  target  echoes. 

Previous  work  was  mainly  concerned  with  large  targets  so  that  the  pulse  lengths  were 
a  fraction  of  the  target  size.  By  contrast  the  targets  investigated  in  this  report  are  small 
anti-personnel  mines  of  sizes  similar  to  the  pulse  length  in  air. 

The  outline  of  this  report  is  as  follows.  Section  2  describes  the  methods  applied 
for  modelling  of  the  GPR  signals.  In  particular,  Section  2.1  presents  the  KT  algorithm 
for  estimating  the  parameters  of  the  exponentially  damped  sinusoid  model.  Section  2.2 
describes  the  two  methods  for  selecting  the  model  order  used  by  the  KT  algorithm.  Section 
3  provides  information  on  the  GPR  system  measurements  and  the  investigated  targets  and 
Section  4  presents  the  results  of  experiments  in  signal  modelling  using  simulated  and  real 
data.  Concluding  remarks  and  recommendations  for  the  future  research  work  on  landmine 
recognition  based  on  the  GPR  echoes  are  presented  in  Section  5. 


2  Theoretical  Background 


The  problem  of  parameter  estimation  of  exponentially  damped  sinusoids  is  very  impor¬ 
tant  in  many  applications.  A  simple  and  widely  used  technique  is  the  Prony  method  [12]. 
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However,  as  it  is  based  on  linear  least-squares,  the  standard  Prony  method  is  highly  sen¬ 
sitive  to  the  additive  measurement  noise.  It  also  requires  knowing  the  number  of  damped 
exponentials  contained  in  the  sum  (i.e.,  the  order  of  the  model)  M  in  advance.  Several 
procedures  to  refine  the  Prony  method  have  been  proposed  [13,  14].  In  particular,  the 
approach  described  in  [13]  based  on  the  principal-component  analysis,  has  been  applied  to 
modelling  the  GPR  backscattered  echoes  (see  [15]).  Unfortunately,  the  author  has  found 
that  the  poles  estimated  using  this  technique  are  highly  unstable. 

The  method  proposed  by  Kumaresan  and  Tufts  [16]  offers  an  improved  performance 
when  applied  to  signals  with  lower  signal-to-noise  ratio  (SNR)  and  compares  favorably 
with  Prony.  It  uses  backward  linear  prediction  and  low-rank  data  matrix  approximation 
based  on  singular  value  decomposition  (SVD).  In  this  report  the  KT  algorithm  is  applied 
to  modelling  the  GPR  echoes  from  our  targets. 

Other  techniques  for  estimating  the  parameters  of  the  assumed  model  from  the  mea¬ 
sured  data  include  ‘‘pencil-of- function”  methods  [17,  18],  maximum  likelihood  approaches 
[19],  techniques  based  on  the  total  least  squares  [20,  21]  and  higher-order  correlations  [22]. 
Barone  at  ai  [23]  proposed  an  approach  for  modelling  time-varying  signals  by  which  the 
observed  interval  is  divided  into  several  short  segments,  each  of  which  is  then  modelled 
separately. 


2.1  Kumaresan  and  Tufts  Algorithm 


Consider  a  data  sequence  x(n)  which  consists  of  M  exponentially  damped  sinusoidal 
signals: 


M 


=  X]  C'fc  exp(5fcn) 
ifc=i 


(1) 


where  Cjt’s  are  non-zero  complex  amplitudes  and  sire  complex  poles  for 

k  —  1, . . .  ,M.  a/c^s  are  the  pole  damping  factors  (a^t  >  0),  and  =  27rfk,  where 
are  the  pole  frequencies.  Let  y(n)  be  an  observed  sequence  corrupted  by  additive  white 
Gaussian  noise  v(n) 


p(n)  =  x(n)  +  v(n)  for  n  =  0, 1,  2, . . . ,  iV  —  1,  N  >  2M. 


For  the  data  sequence  y{n)  backward  linear  prediction  equations  in  matrix  form  are  given 
as: 


y*(i)  y*(2) 

y*(2)  y*(3) 

y*{N-L)  y*{N-L  +  l) 
or,  equivalently, 


y*iL) 

c(l) 

y*(o) 

y*iL  +  i) 

c(2) 

rr  — 

y  +  (1) 

y*{N-i) 

c(L) 

1 

••  1 

* 

AbC  =  -h  (3) 


where  Ab  is  the  {N  -  L)x{L)  backward  data  matrix,  c  is  backward  linear  prediction  filter 
and  denotes  complex  conjugation.  The  above  equation  can  be  written  in  an  augmented 
form  as, 


AU'  =  0 


(4) 
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with  =  (h|Ab)  and  c'  =  where  c'  is  the  prediction-error  filter  and  ‘T” 

denotes  matrix  transpose.  If  the  data  is  noiseless,  the  prediction-error  filter  polynomial 

B{z)  =  1  +  c{l)z~^  +  c(2)z”’^  +  . . .  +  c{L)z~'^  (5) 

has  zeros  at  for  k  =  1,2, provided  that  the  order  of  the  backward  linear 
prediction  equations  in  Eq.(2)  L  is  chosen  such  that  M  <  L  <  N  —  M  [24].  These  M 
zeros  are  the  signal  zeros  and  fall  outside  the  unit  circle  as  the  consequence  of  the  fact 
that  data  is  used  in  the  reverse  direction  to  compute  c'  in  Eq.(2).  If  L  >  M,  polynomial 
B{z)  (5)  has  additional  L  —  M  zeros,  which  are  called  extraneous  zeros.  In  this  case  Eq.(2) 
has  more  than  one  solution  since  the  rank  of  Ab  (or  A{^)  is  M(<  L).  It  has  been  shown 
[24]  that  there  exists  a  unique,  minimum-norm,  solution  to  Eq.(2)  for  which  all  L  —  M 
extraneous  zeros  fall  within  the  unit  circle  provided  that  data  is  noiseless.  This  solution 
corresponds  to  the  least  squares  solution. 

In  this  way  it  is  possible  to  identify  the  M  signal  zeros,  lying  outside  the  unit  circle, 
as  opposed  to  the  L  —  M  extraneous  zeros.  Moderately  large  values  of  L  have  been  found 
essential  in  improving  the  accuracy  of  the  pole  location  estimates  [16].  (Note  that,  for 
L  =  M  and  forward  instead  of  backward  prediction  in  Eq.(2),  this  technique  is  equivalent 
to  the  standard  Prony  method). 

In  the  presence  of  additive  noise  the  least  squares  solution  to  Eq.(2)  is  ill-conditioned 
and  introduces  considerable  perturbations  in  the  solution  vector  c^.  To  alleviate  this 
problem  Kumaresan  and  Tufts  [16]  used  the  truncated  SVD  to  increase  the  SNR  in  the 
data  prior  to  obtaining  the  solution  vector  c.  Namely,  SVD  of  the  data  matrix  Ab  is  given 
as: 

A^  =  U  p  (6) 

where  U  is  an  {N  —  L)x{N  —  L)  dimensional  matrix  of  left  singular  vectors,  V  is  an  LxL 
dimensional  matrix  of  right  singular  vectors  and  S  is  a  diagonal  matrix  of  singular  values 
If  fhe  measurements  are  noiseless,  i.e.,  if  y{n)  =  x(n),  where  x{n)  is  defined  by 
Eq.(l),  then  matrix  Ab  has  rank  M  and  only  M  singular  values  of  Ab  are  nonzero.  In 
the  presence  of  noise  and  model  mismatch  however  the  matrix  Ab  will  have  full  rank,  but 
will  be  “close”  to  a  matrix  Ag  of  the  rank  M,  and  thus  ctm+i  ^  <^m+2  ^  ^  gl  ^ 

0.  Consequently,  Kumaresan  and  Tufts  used  the  optimum  rank  M  approximation  of  A 
obtained  by  setting  L  —  M  smaller  singular  values  of  A  to  zero,  and  computed  the  solution 
vector  c  as, 

M 

c  =  -J2  (7) 

k=l 

where  {cri,(72, . . .  ,crAf}  are  M  largest  singular  values  of  A,  {ui,U2, . . .  ,um}  are  corre¬ 
sponding  left  singular  vectors,  and  {ui,t;2, . . .  are  the  right  singular  vectors. 

Once  signal  poles  are  computed,  complex  amplitudes  Ck  can  be  found  using  linear 
least-squares  approach  which  involves  minimising  the  following  squared-error  sum: 

N  M 

^  H  exp(sfcn)  p  (8) 

n=l  A;=l 

^Least  squares  solution  to  Eq.(2)  can  be  obtained  by  matrix  inversion  using  singular  value  decomposi¬ 
tion. 
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2.2  Choosing  the  Order  of  the  Model 

Equation  Eq.(7)  requires  that  the  order  of  the  model  (or,  simply,  the  model)  M  is 
known.  However,  in  practice  the  correct  model  M  is  not  known  a  priori  and  hcis  to  be 
estimated.  This  estimate  can  be  obtained  using  singular  values  of  the  data  matrix  Au- 
Following  the  discussion  about  the  properties  of  the  singular  values  of  the  matrix  Ab 
corresponding  to  a  noisy  signal  (see  Section  2.1),  a  simple  way  to  adaptively  determine 
the  order  of  the  model  involves  choosing  M  such  that  is  the  first  singular  value 

smaller  than  tai^  where  the  constant  t  is  chosen  from  the  interval  (0, 1)  and  the  singular 
values  are  ordered  in  decreasing  order. 

Several  more  complex  criteria  based  on  the  application  of  information-theoretic  princi¬ 
ples  to  model  selection  are  also  available  (see  [25]  page  548  for  more  details).  The  popular 
methods  are  the  Akaike  information  criterion  (AIC)  [26]  and  the  minimum  description 
length  criterion  (MDL)  [27]  defined,  respectively,  as: 

AIC  =  -^21n(/(y|0))  +  2d  (9) 

MDL  =  -21n(/(y|0))  +  dlnAr.  (10) 

Here  f{y\0)  is  the  maximum  likelihood  function  of  the  observation  vector  y,  6  is  the 
parameter  set  of  the  model  in  Eq.(l)  and  d  is  the  number  of  free  parameters  in  6.  N  in 
Eq.(lO)  denotes  the  total  number  of  observations  in  y.  Both  AIC  and  MDL  select  the 
model  which  minimises  the  corresponding  equation  above. 

In  this  report  we  apply  an  approach  for  determining  AIC  and  MDL  which  is  matched 
to  the  methods  based  on  SVD  and  thus  can  be  used  along  with  the  KT  algorithm  involving 
very  little  extra  computation  [28]. 


3  GPR  Measurements  and  Targets 

In  the  experiments  we  used  an  FR-127-MSCB  MK2  impulse  GPR  system  developed  by 
CSIRO  [29].  The  system  collects  127  echoes,  or  soundings,  per  second  and  each  sounding 
is  composed  of  512  samples  of  12  bit  accuracy.  The  sounding  range  is  variable  between  4 
ns  and  32  ns.  The  system  uses  bistatic  bow-tie  antennas  with  the  transmit  pulse  created 
by  a  fast-recovery  diode  at  the  transmit  antenna.  Centre  frequency  and  bandwidth  of  the 
transmitted  pulse  can  be  changed  by  changing  the  antenna.  In  this  report  we  utilise  the 
soundings  obtained  using  the  antenna  with  1.4  GHz  centre  frequency  and  the  bandwidth 
of  1.136  GHz. 

Six  targets  were  measured  (see  Table  1).  Three  of  these  targets  were  surrogate  land¬ 
mines  developed  by  the  DSTO  countermine  research  project  [30].  These  targets,  denoted 
as  ST-AP(l),  ST-AP(2)  and  ST-AP(3),  are  the  surrogate  anti-personnel  mines  modelled 
after  the  M14,  PMN  and  PMN2  blast  mines,  respectively.  The  PMN  and  PMN2  are  anti¬ 
personnel  mines  with  non-metallic  casings  which  have  been  encountered  by  the  Australian 
troops  attached  to  the  United  Nations  de-mining  operations  in  Cambodia  and  Afghanistan. 
The  M 14  is  a  very  small  anti-personnel  mine  with  almost  no  metal  content  and  very  small 
size.  As  such  it  is  a  very  difficult  target  for  GPR  landmine  detection.  Three  simple  test 
targets,  PVC  (poly-vinyl  chloride)  cylinders,  were  also  measured. 
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List  of  Targets 


Target  Designation 

Target  Description 

Target  Orientation 

PVC  Test  Target  1 

Solid  PVC  Cylinder, 

5.2  cm  in  diameter  and  10.1  cm  high 

Not  critical, 
target  symmetric 

PVC  Test  Target  2 

Solid  PVC  Cylinder, 

10.35  cm  in  diameter  and  5.05  cm  high 

Not  critical, 
target  symmetric 

PVC  Test  Target  3 

Solid  PVC  Cylinder, 

15.35  cm  in  diameter  and  5.05  cm  high 

Not  critical, 
target  symmetric 

ST-AP(l) 

M14  Surrogate  Anti-Personnel  Mine, 
made  from  5.2  cm  diameter  plastic  pipe 
4.2  cm  high,  filled  with  paraffin  wax 
and  some  small  metal  parts. 

Not  critical, 
target  symmetric 

ST-AP(2) 

PMN  Surrogate  Anti-Personnel  Mine, 
made  from  11.8  cm  diameter  plastic  pipe 

5  cm  high,  filled  with  paraffin  wax 
and  some  small  metal  parts 

Trigger  mechanism  of 
the  target  perpendicular 
to  the  direction  of  sweep 

ST-AP(3) 

PMN2  Surrogate  Anti-Personnel  Mine, 
made  from  11.5  cm  plastic  pipe 

5.3  cm  high,  filled  with  paraffin  wax 
and  some  small  metal  parts 

Trigger  mechanism  of 
the  target  perpendicular 
the  direction  to  of  sweep 

Table  1:  Targets  used  in  the  experiments. 


Targets  were  buried  in  dry  sand  approximately  0.5  cm  deep.  The  antenna  was  hung 
from  a  track  along  which  it  could  be  driven  at  constant  velocity  by  a  stepper-motor  at  a 
height  of  approximately  0.5  cm  above  the  surface  of  the  sand  (see  Figure  1).  Each  target 
was  placed  such  that  the  antenna  would  cross  its  centre.  The  antenna  was  oriented  so 
that  the  transmitter  crossed  the  target  first,  followed  by  the  receiver  and,  in  both  cases, 
the  antenna  electric  polarisation  was  perpendicular  to  the  direction  of  travel.  This  results 
in  the  most  complex  target  image  and  highlights  the  internal  structure  of  the  targets. 
The  antenna  size  and  separation  is  of  the  order  of  the  size  of  the  targets.  The  velocity 
of  the  antenna  was  adjusted  to  give  soundings  0.1  cm  apart  and  each  measurement  set 
comprised  approximately  1200  backscattered  GPR  soundings.  These  measurements  were 
initiated  and  terminated  well  off  the  target  so  that  a  large  amount  of  data  was  available 
with  no  target  present.  The  time  range  of  the  GPR  soundings  was  6  ns  which  corresponds 
to  a  sampling  interval  of  At  =  0.0117  ns. 

Shallow-target  GPR  soundings  measured  using  a  bistatic  polar  antenna  are  normally 
obscured  by  the  antenna  cross-coupling  and  surface  return.  These  effects  are  considered 
to  be  approximately  constant  over  a  number  of  soundings  and  may  be  crudely  removed  by 
a  process  called  background  subtraction.  In  this  process,  the  average  of  neighbouring 
echoes  over  a  small  area  with  no  target  present  is  usually  subtracted  from  the  time  domain 
soundings  with  the  target  present.  An  example  of  the  GPR  backscattered  signals  recorded 
in  sand  across  the  surrogate  landmine  ST-AP(2)  is  shown  in  Figure  2.  Figure  2  (a)  shows 
the  raw  measured  signatures.  Figure  2  (b)  shows  the  GPR  signatures  processed  using 
background  removal  with  A/fe  =  1.  The  position  of  the  target  is  clearly  visible  in  this 
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image. 

The  data  files  used  our  in  experiments  are  listed  in  Appendix  B. 


4  Results  and  Discussion 

To  demonstrate  the  performance  of  the  KT  algorithm  we  first  apply  it  to  modelling  a 
simulated  signal  composed  of  two  exponentially  damped  sinusoids  in  noise.  The  complex 
conjugate  poles  of  this  signal  are  Si/2  =  —0.02  ±  j0.18  and  53/4  =  —0.01  ±  jO.l,  and  the 
corresponding  amplitudes  are  C1/2  =  2000  ±  j3  and  C^/4  =  100  ±  j3000.  The  order  of  the 
model  is  M  =  4.  The  signal  is  contaminated  by  additive  zero-mean  white  Gaussian  noise 
with  variance  2a^.  Peak  signal-to- noise  ratio  (SNR)  used  in  this  example  is  set  to  15  dB. 
This  simulated  signal  contaminated  by  noise  is  shown  in  Figure  3. 

The  KT  algorithm  was  applied  to  model  this  noise  contaminated  signal  with  the  order 
of  the  backward  linear  prediction  equations  in  Eq.(2)  set  to  L  =  50.  Two  approaches  de¬ 
scribed  in  Section  2.2  were  used  to  determine  the  order  of  the  model  M.  The  information- 
theoretic  method  [28]  predicted  M  =  6,  whereas  simple  thresholding  of  the  singular  values 
of  the  data  matrix  determined  the  order  as  M  =  38,  where  the  parameter  t  was  set  to 
0.050.  That  is,  in  this  case,  when  the  data  are  truly  derived  from  a  finite  order  autore¬ 
gressive  (AR)  model,  very  good  estimate  of  the  model  order  M  was  obtained  by  using  the 
information-theoretic  criterion,  as  compared  to  the  simple  thresholding  of  the  resulting 
singular  values.  However,  as  will  be  seen  below,  the  information-theoretic  approach  usu¬ 
ally  does  not  perform  as  well  when  the  data  model  is  not  exact,  i.e.,  when  it  is  used  with 
the  real  signals  [25]. 

The  zeros  of  the  prediction  error  filter  polynomial  B{z)  computed  for  the  simulated 
signal  using  the  KT  algorithm  and  the  information-theoretic  criterion  for  the  model  se¬ 
lection  are  shown  in  Figure  4.  The  four  signal  zeros  that  fall  outside  the  unit  circle  are 
clearly  visible,  whereas  the  extraneous  zeros  are  within  the  unit  circle. 

The  approximated  signal  obtained  by  modelling  the  noise  contaminated  data  in  Figure 
3  using  the  KT  algorithm  along  with  the  original  (noiseless)  data  is  shown  in  Figure  5  (a). 
For  comparison.  Figure  5  (b)  shows  the  signal  obtained  by  modelling  the  same  data  using 
the  Prony  method.  The  parameters  of  the  Prony  method  were:  the  order  of  the  model 
M  —  6,  the  spacing  (stride)  =  10  and  no  signal  sub-sampling  has  been  applied  (see  [15] 
for  more  details  on  the  parameters  of  the  Prony  algorithm).  Table  2  presents  the  true  and 
the  estimated  signal  poles  obtained  by  using  the  KT  and  the  Prony  algorithms. 

The  parameters  of  the  Prony  algorithm,  M  and  A:,  used  in  this  example  are  those  for 
which  the  resulting  estimated  poles  and  residues  were  closest  to  their  true  values.  For 
other  combinations  of  M  and  k  the  obtained  estimates  were  much  worse.  The  intention 
was  to  keep  the  parameter  M  as  close  as  possible  to  its  true  value  and  to  obtain,  at  the 
same  time,  good  estimates  of  poles  and  residues.  As  can  be  seen  from  Table  2  for  the 
above  parameters  two  real  poles  and  two  pairs  of  complex-conjugate  poles  were  obtained 
using  the  Prony  method.  We  consider  real  poles  as  spurious  and  not  really  important 
and  treat  the  conjugate-complex  poles  as  our  result.  It  is  interesting  to  note  that  when 
the  order  of  the  model  was  set  to  M  =  4  only  two  complex-conjugate  poles  resulted, 
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Simulation  Results 


True  poles 

Estimated  poles 

KT 

Prony 

-0.02  ±j0.18 
-0.01  ±  jO.lO 

-0.0196  ±  j0.1796 
-0.0079  ±  j0.1007 

0.0036  ±  j0.1946 
-0.098  ±  j0.0916 
0.019 
-0.0016 

True  residues 

Estimated  residues 

KT 

Prony 

2000  ±  j3 

100  ±  j3000 

2158.8  ±  j45.9 
193.6  ±  j2596.5 

100.3  ±  jl37.9 
-879.1  ±  j2722.6 
-8.9 

227.6 

Table  2:  The  true  parameters  of  the  simulated  signal  and  the  parameters  estimated  using 
the  KT  and  the  Prony  algorithms. 


along  with  two  spurious  real  poles.  We  also  mention  that  the  estimates  obtained  by  using 
the  Prony  method  change  randomly  when  the  parameters  M  and  k  are  changing,  and, 
though  it  may  be  possible  to  get  an  estimate  as  good  as  the  one  shown  in  Figure  5  (b), 
the  associated  search  procedure  is  usually  long  and  painstaking.  In  comparison  the  KT 
algorithm  is  much  more  robust  and  outperforms  the  Prony  method  when  applied  to  this 
highly  contaminated  data. 

The  programs  for  modelling  the  GPR  signatures  using  the  KT  algorithm  and  the 
Prony  approach  are  implemented  using  Matlab  Numeric  Computation  and  Visualisation 
Software.  These  programs  are  listed  in  Appendix  C. 

We  now  use  the  KT  algorithm  to  estimate  the  model  parameters  of  the  backscattered 
GPR  echoes  corresponding  to  the  targets  presented  in  Table  1.  Since  the  GPR  [29]  is 
designed  in  such  a  way  that  its  soundings  are  over-sampled,  we  first  sub-sample  these 
signals  by  using  the  appropriate  low-pass  filter  and  decimating  them  by  a  factor  S  =  2 
(5  =  3  gives  very  similar  results).  It  has  been  observed  that  more  exact  pole  estimates 
are  obtained  when  the  signal  is  sub-sampled  in  this  way.  The  reason  is  that  most  of  the 
energy  of  the  originally  acquired  signal  is  concentrated  in  the  lower  part  of  the  spectrum, 
and,  by  filtering  and  decimating  the  signal  in  this  way,  useful  and  interesting  part  of  the 
spectrum  is  distributed  fully  around  the  unit  circle  and  is  not  restricted  to  a  small  angular 
interval.  This  is  desirable,  since  the  resolution  of  the  signal  poles  is  increased  and  their 
more  accurate  estimation  by  the  KT  algorithm  is  enabled. 

Increasingly  better  signal  fits  are  obtained  as  the  order  of  the  model  M  increases.  In 
this  case,  however,  the  estimated  signal  poles  tend  to  be  more  and  more  dispersed.  For 
this  reason  our  intention  was  to  choose  the  minimum  M  for  which  the  signal  is  still  rea¬ 
sonably  well  approximated.  When  applied  to  the  real  data,  the  model  selection  based  on 
thresholding  the  singular  values  of  the  data  matrix  gave  much  better  results  compared  to 
the  information-theoretic  criteria.  Examples  of  modelling  the  GPR  signals  reflected  from 
the  surrogate  landmines  ST-AP(2)  and  ST-AP(3)  using  the  KT  algorithm  and  the  method 
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based  on  thresholding  the  singular  values  of  the  data  matrix  are  shown  in  Figures  6  and  7, 
respectively,  where  X  =  90  (note  that  these  are  raw  GPR  signatures,  i.e.,  without  back¬ 
ground  removal).  The  thresholding  method  resulted  in  M  =  10  and  M  =  12  for  these  two 
targets  for  t  =  0.025.  Prom  Figures  6  and  7  it  can  be  seen  that  the  approximated  signals 
obtained  by  using  the  above  model  orders  were  very  good.  In  contrast,  the  information- 
theoretic  criteria  [28]  gave  M  =  88  for  ST-AP(2)  and  M  =  90  for  ST-AP(3)  which  were 
obviously  too  big.  For  these  reaisons,  in  this  report  we  exclusively  applied  the  thresholding 
method  to  determine  M  for  modelling  the  measured  signatures  of  our  targets. 

The  models  estimated  using  the  thresholding  technique  were  in  most  cases  insensitive 
to  changing  the  parameter  t  within  a  range  of  values.  For  example,  for  the  above  signals, 
where  we  set  t  =  0.025,  this  range  was  from  0.010  to  0.040,  approximately.  This  stability  of 
the  model  estimates  for  varying  values  of  t  within  a  given  range  was  used  as  a  confirmation 
that  the  value  of  the  parameter  {t  =  0.025  in  this  case)  was  appropriately  chosen. 

In  order  to  fully  characterise  each  target  the  signal  poles  are  computed  using  a  series 
of  neighbouring  signatures  recorded  over  each  target  type.  The  poles  estimated  from 
neighbouring  signatures  are  similar  but  not  necessarily  the  same.  The  number  of  poles 
may  also  vary  from  one  signal  to  the  other.  This  spatial  variation  of  signal  poles  can  be 
caused  by  the  internal  structure  of  the  target,  the  characteristics  of  the  imaging  system 
(size  of  the  target  vs.  the  distance  between  the  transmit  and  receive  antennas),  the 
influence  of  noise,  etc. 

The  poles  computed  for  the  targets  in  Table  1  are  projected  in  two-dimensional  a-f 
plane  (see  Figures  8  and  9).  Figure  8  (a)  shows  the  poles  corresponding  to  the  three 
surrogate  landmines  and  Figure  8  (b)  shows  the  poles  computed  for  the  three  PVC  cylin¬ 
ders  where  the  raw  target  signatures  were  used  for  modelling.  The  poles  corresponding  to 
different  target-types  form  clusters  in  the  a-f  space.  It  is  expected  that  the  positions  of 
such  clusters,  which  are  clearly  marked  in  Figures  8  (a)  and  (b),  can  be  used  for  target 
identification.  For  comparison,  Figure  10  shows  the  poles  computed  from  the  (neighbour¬ 
ing)  backscattered  GPR  echoes  corresponding  to  the  case  of  no  target  present  (i.e.,  the 
background  signals). 

Since  the  GPR  signals  are  obscured  by  the  clutter  it  is  also  reasonable  to  apply  the 
background  removal  technique  as  explained  in  Section  3  to  compute  the  target  signatures. 
An  important  question  is  how  to  choose  the  particular  no-target-present  signals  to  be 
subtracted.  In  our  experiments  the  signals  to  be  subtracted  from  the  raw  target  echoes 
are  chosen  randomly  from  the  background  signals  within  the  target-specific  measurement 
sets.  This  process  is  repeated  for  each  target-type  separately.  Figures  9  (a)  and  (b)  show 
the  poles  resulted  from  modelling  such  background  removed  signals  corresponding  to  the 
surrogate  landmines  and  the  PVC  cylinder,  respectively.  It  can  be  seen  that,  though 
the  poles  of  the  background-removed  signals  corresponding  to  different  targets  still  form  a 
number  of  non-overlapping  and  partially  overlapping  clusters,  they  are  also  more  dispersed 
as  compared  to  the  poles  computed  from  the  raw  signatures  (Figures  8  (a)  and  (b)).  This 
effect  is  expected  since  the  background  subtraction  process  effectively  decreases  the  SNR 
of  the  modelled  signals.  For  this  reason  the  parameter  t  in  the  method  for  the  model 
order  estimation  was  set  to  t  =  0.080  for  modelling  the  background-removed  signatures, 
as  opposed  to  the  value  of  t  ^  0.025  which  was  used  when  the  raw  target  signatures  were 
approximated. 
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The  ground  and  target  returns  comprising  the  GPR  signal  can  have  different  time 
shifts  with  respect  to  the  beginning  of  the  mecisurement  (i.e.,  sampling).  This  time  shift 
can  be  adjusted  on  the  GPR  system  at  the  time  when  the  measurements  are  taken,  so  as 
to  take  into  account  the  measurement  conditions,  such  height  of  the  GPR  over  the  ground, 
etc.  If  the  ground  surface  is  relatively  flat  all  signatures  comprising  one  measurement  set 
will  have  similar  time  shifts.  This  shift  may  vary  from  one  measurement  set  to  the  other. 
Since  the  modelling  technique  we  apply  is,  in  general,  not  shift  invariant,  the  estimated 
poles  usually  change  as  the  time  shift  changes.  To  be  able  to  compare  the  poles  extracted 
from  the  signatures  that  belong  to  different  files,  the  modelled  signals  should  be  aligned 
so  as  to  have  the  same  reference  point.  An  investigation  is  needed  into  the  approaches 
for  defining  this  starting  point  automatically  in  a  consistent  manner  for  various  types  of 
signatures. 

To  compute  the  target  poles  in  the  experiments  described  in  this  report  we  used  the 
portion  of  the  original  512-samples-long  GPR  soundings  comprising  the  samples  from 
Nffeg  =  20  to  Nend  ==  400.  Here  the  reference  point  N^eg  was  the  same  for  all  data  files 
since  in  our  experiments  all  measurement  sets  have  been  taken  so  as  to  have  the  same  time 
shift.  Samples  from  N  =  401  to  AT  =  512  have  been  discarded  since  they  do  not  contain 
any  significant  signal  variation. 

The  position  of  the  target  in  a  set  of  measurement  is  determined  by  visual  inspection 
(see  Figure  2).  This  process  can  also  be  automated,  since  the  variances  of  the  GPR 
echoes  reflected  from  the  target  are  substantially  bigger  than  those  corresponding  to  the 
background  signals,  in  particular  when  the  background-subtracted  signatures  are  used. 
Consequently,  a  statistical  test  based  on  the  variance  of  the  signal  might  be  applied  to 
determine  whether  there  is  a  target  present,  and  then  poles  would  be  extracted  from  these 
signals.  It  is  expected  that  the  feasibility  of  this  and  related  approaches  will  be  examined 
in  some  detail  in  the  future. 


5  Conclusions  and  Recommendations 

We  have  investigated  an  approach  to  modelling  the  backscattered  GPR  echoes  from 
various  types  of  buried  non-metallic-cased  surrogate  landmines  and  PVC  cylinders  using 
the  KT  algorithm.  Both  raw  and  background-removed  target  signatures  have  been  con¬ 
sidered.  It  has  been  shown  that  the  poles  extracted  in  this  way  form  clusters  which  are 
characteristic  of  a  given  target-type.  Consequently,  it  is  expected  that  such  poles  can  be 
used  for  the  recognition  of  landmines. 

The  KT  algorithm  has  been  applied  to  synthetic  data  contaminated  by  white  Gaussian 
noise  and  to  experimental  data.  In  modelling  the  synthetic  signal  the  performance  of  the 
KT  method  has  been  compared  to  the  Prony  algorithm  as  applied  in  [15]  with  more  stable 
results. 

Two  diflFerent  approaches  for  selecting  the  order  of  the  model  for  the  KT  algorithm 
have  been  investigated:  an  information-theoretic  and  an  approach  based  on  the  threshold¬ 
ing  of  singular  values  of  the  data  matrix.  While  the  information-theoretic  criterion  was 
successful  for  the  simulated  data  which  were  computed  using  the  the  exact  signal  model, 
the  simpler  thresholding  approach  performed  better  when  applied  to  the  real  data.  The 


9 


DSTO-TR-0680 


poles  computed  using  the  KT  algorithm  and  this  latter  approach  were  stable  for  a  range 
of  threshold  values. 

The  results  presented  in  this  report  correspond  to  only  one  soil  type  (dry  sand)  where 
the  signals  were  measured  using  1.4  GHz  antenna.  It  is  well  known  that  propagation  of 
electromagnetic  waves  depends  on  the  properties  of  the  medium,  and  we  can  expect  to 
obtain  different  results  for  the  soil  types  with  varying  dielectric  constants.  For  example, 
higher  dielectric  constants  introduce  higher  losses,  so  the  returned  signal  is  weaker.  The 
dielectric  constants  also  increase  with  the  moisture  content  of  the  soil.  It  remains  to 
be  seen  how  different  conditions,  i.e.,  soil  types,  depths  at  which  the  target  is  buried, 
different  target  aspects,  etc.,  affect  the  pole  estimates,  and,  what  antenna  bandwidths 
provide  most  information  for  the  recognition.  The  use  of  balanced  bridge  antennas  should 
also  be  considered. 

Further  research  is  needed  to  determine  whether  or  not  to  use  the  background-removed 
signatures,  and  how  to  define  the  background  signal  to  be  subtracted.  A  viable  approach 
to  this  problem  would  be  to  use  the  assumption  that,  when  the  raw  target  signatures 
are  modelled,  some  of  the  resulting  poles  correspond  to  the  antenna  cross-coupling  and 
ground  return,  whereas  other  poles  represent  the  part  of  the  signal  returned  from  our 
target.  Therefore,  the  background  signal  might  be  synthetised  for  each  signature  using 
the  poles  corresponding  to  the  ground  return.  These  synthetised  background  signals  would 
then  be  subtracted  from  the  raw  signals.  In  this  way  some  of  the  poles  originally  obscured 
by  the  clutter  could  be  extracted.  It  is  expected  that  detection/recognition  of  small  targets 
could  be  improved  in  this  way. 

We  are  also  looking  into  developing  a  pattern  recognition  system  for  buried  targets 
based  on  the  pole  estimates  and  we  will  be  fusing  this  information  with  the  signals  obtained 
from  a  metal  detector. 

Following  is  the  summary  of  our  recommendations  for  the  future  work  on  detec¬ 
tion/recognition  of  landmines: 

•  It  is  necessary  to  investigate  target  signatures  for  various  soil  types  and  conditions 
to  determine  the  stability  of  poles. 

•  Signals  recorded  using  different  centre  frequency  antennas  need  to  be  examined  to 
determine  the  frequency  range  optimal  for  detection  and  identification  of  landmines. 

•  It  is  important  to  undertake  further  research  on  choosing  the  most  suitable  back¬ 
ground  removal  method.  In  addition,  approaches  to  determining,  in  a  consistent 
manner,  the  part  of  the  original  512-samples-long  GPR  soundings  to  be  used  in 
the  KT  modelling  for  a  range  of  targets  and  environment  conditions,  need  to  be 
investigated. 

•  Techniques  for  target  detection  using  the  signal  variance  and  other  approaches  need 
to  be  investigated  in  order  to  locate  the  target. 

•  Since  the  GPR  has  difficulty  discriminating  between  the  small  targets  and  clutter,  it 
is  envisaged  that  much  better  identification  results  could  be  obtained  by  fusing  the 
information  obtained  from  the  EMI  metal  detector  and  that  from  the  GPR  alone. 
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In  parallel  with  the  above  it  is  necessary  to  define  the  most  appropriate  pattern 
recognition  method  for  landmine  identification  based  on  the  estimated  poles. 


References 

1.  D.  J.  Daniels,  D.  J.  Gunton,  and  H.  F.  Scott.  Introduction  to  subsurface  radar.  Proc. 
lEE,  135  (F)(4):278-317,  1988. 

2.  D.  J.  Daniels.  Surface  penetrating  radar.  Radar,  sonar,  navigation  and  avionics  series 
6.  The  Institution  of  Engineers,  London,  UK,  1996. 

3.  C.  E.  Baum,  E.  J.  Rothwell,  K.-M.  Chen,  and  D.  P.  Nyquist.  The  singularity  expansion 
method  and  its  appliaction  to  target  id  entification.  Proc.  IEEE,  79(10):1481-1491, 
October  1991. 

4.  L.  C.  Chen,  D.  L.  Moffatt,  and  L.  Peters  Jn.  A  characterization  of  subsurface  radar 
targets.  Proc.  IEEE,  67(7):951-1000,  July  1979. 

5.  L.  C.  Chan,  L.  Peters  Jr.,  and  D.  L.  Moffatt.  Improved  performance  off  a  subsurface 
radar  target  identification  system  through  antenna  design.  IEEE  Trans.  Antenna 
Propag.,  29(2):307-311,  March  1981. 

6.  E.  M.  Kennaugh.  The  K-pulse  concept.  IEEE  Trans.  Antennas  Propag.,  29(2);327- 
331,  March  1981. 

7.  K.-M.  Chen,  D.  P.  Nyquist,  E.  J.  Rothwell,  L.  L.  Webb,  and  B.  Drachman.  Radar 
target  discrimination  by  convolution  of  radar  return  with  extinction-pulses  and  single¬ 
mode  extraction  signals.  IEEE  Trans.  Antennas  Propag.,  34(7):896-904,  July  1986. 

8.  E.  J.  Rothwell,  K.-M.  Chen,  and  D.  P.  Nyquist.  Extraction  of  natural  frequencies  of  a 
radar  target  from  a  measured  response  using  E-pulse.  IEEE  Trans.  Antennas  Propag., 
35(6):715-720,  June  1987. 

9.  Y.  Bissessur  and  R.  N.  G.  Naguib.  Buried  plant  detection:  A  Volterra  series  modelling 
approach  using  atrificial  neural  networks.  Neural  Networks,  9(6):1045-1060,  1995. 

10.  J.  W.  Brooks  and  M.  W.  Maier.  Applicatiion  of  system  identification  and  neural 
networks  to  classification  of  land  mines.  In  Proc.  EUREL  Int.  Conf.,  7-9  Oct.  1996, 
EICC,  Edinburgh,  UK,  pages  46-50,  1996. 

11.  H.  C.  Strifes,  K.-M.  Chen,  S.  Abramson,  B.  Brusmark,  and  G.  C.  Gaunaurd.  Signature 
features  in  time-frequency  of  simple  targets  ectracted  by  ground  penetrating  radar.  In 
Proc.  IEEE  Intrnational  Radar  Conf.,  8-11  May,  Virginia,  1995. 

12.  M.  L.  Blaricum  and  R.  Mitra.  A  technique  for  extracting  the  poles  and  residues 
of  a  system  directly  from  its  transient  response.  IEEE  Trans.  Antennas  Propag., 
23(6):777-781,  November  1975. 

13.  M.  L.  Blaricum  and  R.  Mitra.  Problems  and  solutions  associated  with  Prony’s  method 
for  processing  transient  data.  IEEE  Trans.  Antennas  Propag.,  26(1):174-182,  January 
1978. 


11 


DSTO-TR-0680 


14.  R.  Kumaresan,  D.  W.  Tufts,  and  L.  L.  Scharf.  A  Prony  method  for  noisy  data: 
chosing  the  signal  components  and  selecting  the  order  in  exponential  signal  models. 
Proc.  IEEE,  72(2):230-233,  1984. 

15.  D.  Carevic,  M.  Craig,  and  I.  Chant.  Modelling  GPR  echoes  from  land  mines  using 
linear  combinations  of  exponentially  damped  sinusoids.  In  Proc.  SPIE  AeroSense  ^97, 
Orlandoj  USA,  April  1997. 

16.  R.  Kumaresan  and  D.  W.  Tufts.  Estimating  the  parameters  of  exponentially  damped 
sinusiods  and  pole-zero  modeling  in  noise.  IEEE  Trans.  Acoustic,  Speech,  Signal  Proc., 
30(6):833-840,  1982. 

17.  Y.  Hua  and  T.  K.  Sarkar.  Generalized  pencil-of-function  methods  for  extracting  poles 

of  an  em  system  from  its  transit  reponse.  IEEE  Trans.  Antennas  Propag.,  37(2)  :229- 
234,  February  1989.  - 

18.  Y.  Hua  and  T.  K.  Sarkar.  Matrix  pencil  metod  for  estimating  parameters  of  expo¬ 
nentially  damped/undamped  sinusiods  in  noise.  IEEE  Trans.  Acoustic,  Speech,  Signal 
Proc.,  38(5):814-824,  May  1990. 

19.  Y.  Bresler  and  A.  Macovski.  Exact  maximum  likelihood  parameter  estimation  of 
superimposed  exponentioal  signals  in  noise.  IEEE  Trans.  Acoustic,  Speech,  Signal 
Proc.,  34(5):  1081-1089,  October  1986. 

20.  M.  A.  Rahman  and  K.-B.  Yu.  Total  least  squares  approach  for  frequency  estimation 
using  linear  prediction.  IEEE  Trans.  Acoustic,  Speech,  Signal  Proc.,  35(10):1440-1454, 
October  1987. 

21.  R.  Carriere  and  R.  L.  Moses.  High  resolution  radar  target  modeling  using  a  modified 
Prony  method.  IEEE  Trans.  Antennas  Propag.,  40(1):  13-18,  January  1992. 

22.  D.  P.  Ruiz,  M.  A.  Carrion,  A.  Gallego,  and  A.  Medouri.  Parameter  estimation  of 
exponentially  damped  sinusiods  using  a  higher  order  correlation-based  approach.  IEEE 
Trans.  Signal  Proc.,  43(ll):2665-2667,  November  1995. 

23.  P.  Barone,  E.  Massaro,  and  A.  Polichetti.  The  segmented  Prony  method  for  the 
analysis  of  non-stationary  time  series.  Astron.  Astrophys.,  209:435-444,  1989. 

24.  R.  Kumaresan.  On  the  zeros  of  the  linear  prediction-error  filter  for  deterministic 
signals.  IEEE  Trans.  Acoustic,  Speech,  Signal  Proc.,  31(l):217-220,  1983. 

25.  C.  W.  Therrien.  Discrete  random  signals  and  statistical  signal  processing.  Prentice 
Hall,  Englewood  Cliffs,  NJ  07632,  1992. 

26.  H.  Akaike.  A  new  look  at  the  statistical  identification.  IEEE  Trans.  Automat.  Contr., 
19:716-723,  December  1974. 

27.  J.  Rissanen.  Modeling  by  shortest  data  description.  Automatica,  14:465-471,  1978. 

28.  V.  U.  Reddy  and  L.  S.  Biradar.  SVD-based  information  theoretic  criteria  for  detection 
of  the  number  of  damped/undamped  sinusoids  and  their  performance  analysis.  IEEE 
Trans.  Signal  Proc.,  41(9):2872-2881,  1993. 


12 


DSTO-TR-0680 


29.  W.  Murray,  Williams  C,  J,  and  J.  T.  A.  Pollock.  A  high  resolution  radar  for  mine 
detection.  In  Proc.  EUREL  Int.  Conf.,  EICC,  Edinburgh,  UK,  7-9  Oct.  1996. 

30.  B.  B.  Y.  Wong,  I.  Chant,  G.  N.  Crisp,  K.  Kappra,  K.  Strugess,  A.  Rye,  and  K.  Sher- 
bondy.  Suggested  soil  characterisation  techniques  and  surrogate  targets  for  ultra-wide¬ 
band  radar  mine  detection  experiments.  In  Proc.  SPIE  AeroSense  ’97,  Orlando,  USA, 
April  1997. 


13 


DSTO-TR-0680 


Appendix  A 
Figures 


Figure  1:  The  GPR  experiment  in  progress.  Pictured  is  the  antenna  track  sitting  across 
the  dry  sand  box  with  the  GPR  attached. 
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Figure  2:  GPR  backscattered  signals  taken  in  sand  across  the  surrogate  landmine  target 
ST’‘AP(2),  (a)  original  (raw-data)  soundings;  (b)  background  signal  subtracted  from  the 
data.  The  position  of  the  target  in  the  image  is  indicated. 
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Figure  4:  Zeros  of  the  prediction- error  filter  polynomial  B{z).  The  signal  poles  are  outside 
the  unit  circle  whereas  the  extraneous  zeros  falls  within  the  unit  circle. 
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(a) 


(b) 

Figure  6:  GPR  signature  of  the  surrogate  land  mine  ST-AP(2):  (a)  time-domain  represen¬ 
tation  (normalised  to  the  values  between  1  and  -1);  (b)  frequency- domain  representation. 
Solid  line  represents  the  original  signature  and  dashed  line  corresponds  to  the  signal  mod¬ 
elled  using  the  KT  algorithm.  The  method  based  on  the  thresholding  of  the  singular  values 
of  data  matrix  was  applied  to  determine  the  model  order  with  t  =  0.025. 
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(a) 


(b) 

Figure  7:  GPR  signature  of  the  surrogate  land  mine  ST-AP(3):  (a)  time-domain  represen¬ 
tation  (normalised  to  the  values  between  1  and  -1);  (b)  frequency-domain  representation. 
Solid  line  represents  the  original  signature  and  dashed  line  corresponds  to  the  signal  mod¬ 
elled  using  the  KT  algorithm.  The  method  based  on  the  thresholding  of  the  singular  values 
of  data  matrix  was  applied  to  determine  the  model  order  with  t  =  0.025. 
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(a) 


(b) 

Figure  8:  Poles  of  the  raw  GPR  signatures  corresponding  to  the  test  targets  listed  in  Table 
1  estimated  using  the  KT  algorithm  with  L  =  90;  (a)  surrogate  land  mines:  ‘o^  -  ST- 
AP(1),  -  ST-AP(2),  -  ST-AP(3);  (b)  PVC  cylinders:  ^o’  -  PVC  Test  Target  1, 

-  PVC  Test  Target  2,  -  PVC  Test  Target  3,  The  method  based  on  the  thresholding  of 

the  singular  values  of  data  matrix  was  applied  to  determine  the  model  order  with  t  =  0.025. 
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Damping  coefficients  x  10e+9 


(a) 


Damping  coefficients  x  1 0e+9 


(b) 

Figure  9:  Poles  of  the  background  removed  GPR  signatures  corresponding  to  the  test  tar¬ 
gets  listed  in  Table  1  estimated  using  the  KT  algorithm  with  L  =  90;  (a)  surrogate  land 
mines  background  subtracted:  ^o^  -  ST-AP(l),  ^  -  ST-AP(2),  -  ST-AP(3);  (b)  PVC 

cylinders:  ^o^  -  PVC  Test  Target  1,  -  PVC  Test  Target  2,  -  PVC  Test  Target  3. 

The  method  based  on  the  thresholding  of  the  singular  values  of  data  matrix  was  applied  to 
22  determine  the  model  order  with  t  =  0.080. 


DSTO-TR--0680 


Damping  coefficients  x  1 0e+9 


Figure  10:  Poles  of  the  GPR  backs cattered  signals  corresponding  to  the  case  of  no  target 
present  (background  signals)  estimated  using  the  KT  algorithm  with  L  =  90,  The  method 
based  on  the  thresholding  of  the  singular  values  of  data  matrix  was  applied  to  determine 
the  model  order  with  t  =  0.025. 


DSTO-TR-0680 


Appendix  B 

GPR  Measurement  Sets 


Data  files  containing  the  measurements  used  in  the  experiments  were  obtained  by  FR~127- 
MSCB  MK2  impulse  GPR  system.  (Reference  [29]  provides  detailed  description  of  this 
GPR  system.)  The  1.4  GHz  centre  frequency  antenna  connected  in  non-differential  mode 
was  used  and  the  soil  was  dry  sand.  The  measurement  procedure  and  the  investigated 
targets  are  described  in  Section  3  of  this  report. 

Table  3  shows  the  data  files  and  the  corresponding  targets. 


Target 

File  Name 

ST-AP(l) 

ST-AP(2) 

ST-AP(3) 

/ANT141  /SAND /SURRl /FILEl 
/ANT141/SAND/SURR3/FILE1 
/ANT141/SAND/SURR2/FILE1 

PVC  Test  Target  1 
PVC  Test  Target  2 
PVC  Test  Target  3 

/ANT141/SAND/PVCCYL1/FILE1 
/  ANT141  /SAND /P  VCCYL2/FILE1 
/ANT141/SAND/PVCCYL3/FILE1 

Table  3:  GPR  data  files  used  in  the  experiments. 
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Appendix  C 
Program  Listings 


7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


Function :  [P , R , A] =kt alg (D,L,step,sc,fb, svdt ) 

Models  the  signal  as  a  sum  of  exponetially  damped  sinusois  using  the 
KT  algorithm. 


Input:  D  -  real  matrix  comprising  m  GPR  traces  of  length  n 

(averaged  to  obtain  a  single  trace) 

L  -  maximum  reqired  number  of  poles 

sc  -  if  sc>l  data  is  low-pass  filtered  and  correspondingly 
decimated,  sc  values:  2,  3,  4,  etc. 
fb  1  -  forward  prediction 
-1  -  backward  prediction 
step  -  required  spacing  between  successive  samples 
svdt  -  parameter  for  thresholding  the  singulau:  values 
of  the  data  matrix 

Output:  P  -  complex  vector  of  dimension  dim=30  containing  signal 

poles  and  the  following  data:  deg  (dim),  ste  (dim-1), 
sc  (dim-2)  and  fb  (dim-3) 

R  -  complex  vector  containing  corresponding  residues 
A  -  real  vector  containing  approximation  of  the  input  trace 


Author:  Dragana  Carevic 

DSTO  Tactical  Surveillance  Systems  Division 
tel:  8259  6804 

email :  dragana . carevicOdsto . defence . gov . au 
Date:  6/4/1997 


Kumauresan  and  Tufts  algorithm  explained  in: 

R.  Kumairesan  and  D.  W.  Tufts 

"  Estimating  the  parameters  of  exponentially  damped  sinusiods 
and  pole-zero  modeling  in  noise,  "  IEEE  Trans.  Acoustic, 
Speech,  Signal  Proc.  pp.  833-840,  1982. 

SVD  based  information-theoretic  criteria  for  order  selection 
as  per: 
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y,  V.  U.  Reddy  and  L.  S.  Biradar 

7,  "  SVD-based  information  theoretic  criteria  for  detection 

7,  of  the  number  of  dumped/\mdumped  sinusiods  and  their 

7.  performance  analysis,"  IEEE  Trans  Sig.  Process.  Vol.  41 

7,  No.  9,  September  1993. 

7. 

7,  Uses:  reconst. m 

7. 


function  [P,R,A]  =  ktalg(D,L,step,sc,fb,svdt) ; 


[n,m]=size(D) ; 

[MM,beg3=max(D) ; 

if  m>l 

f  p=mean  (D  (  :  ,  1 :  m)  O  ' ; 
m=l ; 
else 
fp=D; 
end; 

7.  Low-pass  filtering  and  decimation  if  sc>l 
if  sOl 

lh=50;  y.length  of  the  low-pass  fir  filter 
f=decimatf  (fp,sc,lh,  ^firO ; 

[n,l]=size(f) ; 
else 
f=fp; 
end; 

clear  fp; 

ave=mean(f ) ; 
f=f-ave; 

7.  Start  here 
dim=30 ; 

P=zeros(dim,l) ; 

R=zeros(dim,l) ; 

A=zeros(n,l) ; 
nl=n-step’*'L 

7.  Compute  data  matrix  FP 

if  fb==l  7.F0RWARD  PREDICTION 

f print f ( ^  Forward\n ’ ) ; 
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for  k=l : (nl) 

t=k:step: ((L-l)*step+k) ; 

FP(k.:)=f(t)’; 
y(k)=f (L*step+k) ; 
end; 

elseif  fb==-l  ‘/.BACKWARD  PREDICTION 
fprintf ( ’Backward\n’ ) ; 
for  k=l:nl 

t=k+step : step : L*step+k ; 

FP(k,:)=f(t)’; 

y(k)=f(k); 

end; 

end; 

*/,  Compute  singular  value  decomposition  of  FP 
[U,S,V]=svd(FP); 

S(l:L/2.1:L/2); 

J=diag(S) ; 


y.  Apply  SVD-based  information-theoretic  criteria  for 
*/,  detecting  the  number  of  poles 

h=U’*y’;  HN=norm(h,2)“2; 
r=nl; 

for  k=l:L/2 

HN=HN- (h(2*k-l) ) “2- (h(2*k) ) ~2 ; 
sum=r*(l+log(2*pi)+log(HN/r) ) ; 

AIC(k , l)=sum+2* (2*k+l) ; 

MDL(k, l)=sum+(2*k+l)*log(r) ; 
end; 

[mrnin , k] =min ( [AIC , MDL] ) ; 

Cm,kl]=min(mmin) ;  Ml=2*k(kl); 

fprintf  (’Information-theoretic  model  order  is  */,d\n’,Ml); 

%  A  heuristic  criterion  for  determining  model  order 
for  k=3:L 

if  S(k,k)<svdt*S(l ,  1)  */,  0.025  for  sc=2  and  raw  signal  ant  1.4  GHz 

*/,  0.050  for  sc=2  and  backgroimd  removed  data  1.4  GHz 
*/,  0.050  fo  sc=2,  raw  data  3GHz 

break; 

end; 

end; 

if  rem(k,2)==0 
M2=k; 
else 
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M2=k-1; 

end; 

fprintf ( ’Heuristic  model  order  is  7,d\n’,M2); 

*/,  Compute  vector  b 
M=M2; 

S=diag(S) ; 

S=diag(ones(M,l) ./S(l :M)) ; 
b=zeros(L,l) ; 

b=-V(; ,1:M)*S*U(: ,1:M) ’*y’ ; 

PC=zeros(L+l,l) ; 

PC(1,1)=1;  PC(2:L+l,l)=b; 

*/,  Solve  characteristic  equation 
RTS=roots((PC)) ; 

7,  Plot  roots  on  the  unit  circle 
7.f  igure ;  plot  (RTS  grid ; 

7thold  on;  7.  plot  \uiit  circle 
7tW=0 : 2’t‘pi/lOO :  2’*‘pi ; 

7ta=cos(w)+j*sin(w)  ; 

7tplot(imag(a)  ,real(a))  ; 
7.axis([-1.5,l-5,-1.5,l-5])  ; 

7.  Choose  poles  which  fall  outside  unit  circle 

ist=l/step; 

if  fb==-l 

i=find(abs(RTS)>l) ; 

[Kl,K2]=size(i) ; 

K=K1;  if  Kl==l 

K=K2  7tthis  is  final  number  of  poles 
end; 

for  k=l:K 

P(k,l)=-ist*log(RTS(i(k))) ; 
if  abs(imag(P(k)))*step==pi 
P(k,l)=real(P(k,l)); 

end; 

end; 

else 

P(l:L,l)=ist.*log(RTS);K=L; 

end; 

RTS=RTS(i); clear  i; 

[Y,i]=sort(abs(imag(P(l :K,1)))) ; 
PS=zeros(dim,l) ; 
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PS(l:K.l)=(P(i)); 

P(diin,l)=K; 

P(dim-1 ,l)=step; 

P (dim-2, l)=sc; 

P (dim-3, l)=fb; 

PS(dim-3:dim,l)=P(dim-3:dim,l) ; 

P=PS; 

deg=K; 


*/.  Compute  the  residuals  and  signal  approximation 
[res,a]=reconst(f ,P) ; 


R(l:deg)=res;  ‘/.residuals 
A=a+ave ; 

*/,  aproximate  signal 
f o=f+ave ; 

7,  input  signal 

7.  Plot  the  reconstructed  signal  in  time  and  frequency  domains 

TB=6;  7.  time  base  in  nanoseconds  6  for  1.4  GHz  antenna 
7.  10  for  3.0  GHZ  antenna 
TP30=sc*TB/512;  7.in  ns 
FP30=1/TP30;  7.in  GHZ 

tt=(0:TP30: (n-l)*TP30) ;mmax=meoc(fo) ; 

f o=f o/mmax ; Al=A/mmax ; 

figure;  plot(tt,fo, ’b-’) ;  grid; 

hold  on;  plot(tt,Al, ’r — ’);hold  off; 

xlabeK’Time  (ns)’);  ylabel( ’Relative  Amplitude’); 


four(: ,l)=10*logl0(fftshift(abs(fft(fo)))) ; 

four ( : , 2) =10*logl0 (f f tshif t (abs (f f t (Al) ) ) ) ; 

ff=(-FP30/2:FP30/n: (n-l)*FP30/(2*n)); 

figure ;plot(ff,f our (: ,1) ,  ’b-’) ;grid; 

hold  on;  plot(ff  ,fo\ir(:  ,2)  , ’r — ’);  hold  off; 

xlabel  (’Frequency  (GHz)’);  ylabeK ’Magnitude  (dB)’); 

end;  7.  function 
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y^***:^**!ti*:lf^f4iilfi1i*1>i**i>i:^:^^c*:t‘***************************************************** 

7. 

7.  Function:  CP,R,A]=pronalg(D,L,step,sc) 

7. 

7.  Models  the  signal  as  a  sum  of  exponetially  damped  sinusoids 

7.  using  Prony  method.  Prony  signal  modelling  involves  computing 

7.  poles  using  principle  component  analysis. 

7. 


7. 

7.  Input : 

7. 

7. 

7. 

7. 

7. 

7  Output : 

7. 

7. 

7. 

7. 

7. 

7.  Author : 

7. 

7. 

7. 

%  Date : 

7, 


D  -  real  matrix  comprising  m  GPR  traces  of  length  n 
(averaged  to  obatin  a  single  trace) 

L  -  reqired  number  of  poles 

sc  -  if  sc>l  data  is  low-pass  filtered  and  correspondingly 
decimated 

step  -  required  spacing  between  successive  samples 
P  -  complex  vector  of  dimension  dim=30,  containing  signal 
poles  and  the  following  data:  L  (dim),  step  (dim-1) 
auid  sc  (dim-2) 

R  -  complex  vector  containing  corresponding  residues 
A  -  real  vector  containing  approximation  the  input  trace 

Dragana  Care vie 

DSTO  Tactical  Surveillance  Systems  Division 
tel:  8259  6804 

email :  dragana . carevicOdsto . defence . gov . au 
6/4/1997 


7.  Description:  Input  signal  undergoes  ’pole-zero’  analysis  and  is 
7.  approximated  using  a  constant  offset  term  and  several 

7.  damped  sinusoidal  oscillations,  each  treated  as  a  pair 

7.  of  conjugate  complex  exponetials. 

7. 

7.  On  output  poles  x+iy  and  residues  u+iv  describe  a  conjugate 

7.  pair  of  model  terms 

7. 

%  (u+iv) exp [ (x+iy) t]  and  (u-iv)exp[(x-iy)t] . 

7. 

Vt  Predicted  curves  according  to  this  model  are  retiirned  in 

7.  a  real  matrix  A. 

7. 


function  [P,R,A]  =  pronalg(D,L,step,sc) ; 


[n,m]=size(D) ; 
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[MM,beg]=max(D) ; 
if  m>l 

f p=mean (D(:,l:m)’)’; 
m=l ; 
else 
fp=D; 
end; 

y.  Low-pass  filtering  and  decimation  if  sc>l 
if  sc>l 

lh=50;  yiength  of  the  low-pass  fir  filter 
f=decimatf  (fp,sc,lh,  ^firO ; 

Cn,l]=size(f ) ; 
else 
f=fp; 
end; 

clear  fp; 

ave=mean(f  )  ; 
frsf^ave; 

%  Start  here 
dim=30 ; 

P=zeros(dim, 1) ; 

R=zeros(dim,l) ; 

A=zeros(n, 1) ; 
nl=n“Step*L 


%  Eigenvector  solution 
for  k=l:(nl) 

t=k:step:  ((L-D+step+k)  ; 
FP(k,:)=f(t)^ 
y(k)=f (L^step+k) ; 
end; 

FPP(:,l)=y^  FPP(:  ,2:L+1)=FP; 

cov=FPP  ^ ♦FPP . /nl ; 
[VEC,D]=eig(cov); 
incom=zeros(L+l,l) ; 
for  k=l:L+l 
incom(k, l)=D(k,k) ; 
end; 

[1 , 1]  =sort  (incom)  ; 
I=reverse(I) ; 
PC=VEC(:,I(L+1)); 
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7,  Solve  characteristic  equation 
RTS=roots((PC)); 

7.  Plot  root  and  the  unit  circle 
7.  f  igure;  plot  (RTS, ’k*’)  ;grid; 

7.  hold  on;  7.  plot  unit  circle 
7.  w=0:2*pi/100:2*pi; 

7t  a=cos(w)+j*sin(w) ; 

7,  plot(imag(a)  ,real(a) , ’k’)  ; 

7.  axis([-1.5,1.5,-1.5,1.5]); 

7  Compute  signal  poles 

ist=l/step; 

for  k=l:L 

if  imag(RTS(k))==0 

P (k , l)=sign(RTS (k) ) ♦ist*log(abs (RTS (k) ) ) ; 
else 

P(k,l)=ist.*log(RTS(k)); 

end; 

end; 

K=L; 

[Y,i]=sort(abs(imag(P(l:K,l)))) ; 

PS=zeros(dim, 1) ; 

PS(l:K,l)=(P(i)); 

PS(dim,l)=K; 

PS (dim- 1 , 1 ) =step ; 

PS (dim-2, l)=sc; 

P=PS; 


7,  Compute  residuals  and  signal  approximation 

[res,a]=reconst(f ,P) ; 

deg=K; 

R(l:deg)=res;  7.  residuals 
A=a+ave ; 

7,  approximated  signal 

f o=f +ave ; 

7.  input  signal 

7.  Plot  the  reconstructed  signal  in  time  and  frequency  domains 
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A=a+ave ; 
f o=f +ave ; 

TB=6;  */,  time  base  in  nano  seconds  6  for  1.4  GHz  antenna 
7,  10  for  3.0  GHZ  antenna 

TP30=sc*TB/512;  y,in  ns 
FP30=1/TP30;  V.in  GHZ 

tt=(0:TP30: (n-l)*TP30) ;mmax=max(f o) ; 

f o=f o/mmax ; Al=A/mmax ; 

figure;  plot(tt,fo, ’b-’) ;  grid; 

hold  on;  plot(tt,Al, ’r — ’);hold  off; 

xlabeK’Time  (ns)’);  ylabeK ’Relative  Amplitude’); 

four( : , l)=10*logl0(fftshift(abs(fft (fo)))) ; 
four( : ,2)=10*logl0(fftshift(abs(fft(Al)))) ; 
f f = (-FP30/2 : FP30/n ; (n-1 ) ♦FPBO/ (2*n) ) ; 
figure ;plot(ff .four ( : ,1) ,  ’b-’) ;grid; 
hold  on;  plot (ff .four ( : ,2) , ’r — ’);  hold  off; 
xlabel  (’Frequency  (GHz)’);  ylabeK ’Magnitude  (dB)’); 


end;  7.  fimction 
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•/. 

•/. 

•/. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


4c  4c  4c  :|c  :4c  4c  ^  :(c  ]4c  1 ^  4c  ♦  4c  %%  4c  :tc  :((:4c 

Function:  [res,a3=reconst(f ,pol) 

Uses  the  computed  poles  to  compute  the  residues  and  the 
model-based  approximation  of  the  signal. 

Input:  f  -  input  signal  of  length  n 

P  -  complex  vector  of  dimension  dim=60,  containing  signal 
poles  and  the  following  data: 
deg  (i=dim) ,  h  (i=dim-l) . 

Output:  res  -  complex  vector  containing  corresponding  resid 

a  -  real  vector  containing  approximation  of  the  input 
signal  f 

Author:  Dragana  Cairevic 

DSTO  Tactical  Surveillance  Systems  Division 
tel:  8259  6804 

email :  dragana . carevicQdsto . defence . gov . au 
Date:  20/4/1997 

4c  4c  4c  4t  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4(  4c  4c  4t  4c  4c  ♦♦  4c  4c  4c  4c  4c  ♦  4c  4c  4c  4c  4c  ♦♦  4c  4c  4c  4c  4c  4c  4c  4c  4c  *  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  *  4c  4c  4c  4c  4c  4c  4c  4(  4c  4c  4c  4c  4c  4c  4c 


function  [res, a]  =  reconstCf ,pol) ; 


[d,m]=size(pol) ;  %  d  is  the  number  of  poles,  m=l 
deg=pol(d) ; 
step=pol(d-l) ; 

[n,l]=size(f) ; 
a-zeros(n,l) ; 

C=zeros(n,deg) ; 
for  k=l:n 
1=1;  dt=l; 
while (l<=deg) 

r=real(pol(l)) ;  w=imag(pol(l)) ; 
C(k,l)=exp(r4ck4cdt)*cos(w4ck4cdt) ; 

1=1+1; 
if  (w'"=0) 

C  (k ,  1)  =-exp  (r 4ck4cdt )  4c  s  in  ( w4ck*dt )  ; 

1=1+1; 
end;  %\± 
end;  %  while 
end;  %  k  loop 


35 


DSTO-TR--0680 


X=inv(C'*C)^C'*f ;  X  LSE  solution  of  CX=f 
clear  C; 

res=zeros(deg,l) ; 

1-1; 

while (l<=deg) 

if  imag(pol(l))==0  ‘/.real  pole 
res(l)=X(l); 

1=1+1; 

else 

res(l)=X(l)  ./2+X(l+l)  ./2*j  ;  7,  conjugate  complex  amplitudes 
res(l+l)=X(l) ./2-X(l+l) ./2*j ; 

1=1+2; 
end;  y,if 
end;  7,  while 
clear  X; 

7.  Compute  the  reconstructed  signal 
a=zeros(n,l) ; 
for  k=l :n 
for  1=1: deg 
if  abs(res(l)'"=0) 

a(k)=a(k)+res(l)*exp(pol(l) *k) ; 
end; 

end;  7.  1  loop 
end;  7.  k  loop 

end;  7«f  unction 
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